function printReconstructionContent(modelFolder,propertiesFolder,reconVersion,numWorkers)
% This function creates text files containing all reactions and metabolites
% in the reconstruction resource.
%
% USAGE
%   printReconstructionContent(modelFolder,propertiesFolder,reconVersion,numWorkers)
%
% INPUTS
% modelFolder         Folder with reconstructions to be printed
% propertiesFolder    Folder where the computed stochiometric and flux
%                     consistencies will be stored
% reconVersion        Name assigned to the reconstruction resource
% numWorkers          Number of workers in parallel pool
%
%   - AUTHOR
%   Almut Heinken, 07/2020

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end
solver = CBT_LP_SOLVER;

if numWorkers > 0
    % with parallelization
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        parpool(numWorkers)
    end
end
environment = getEnvironment();

dInfo = dir(modelFolder);
models={dInfo.name};
models=models';
models(~(contains(models(:,1),{'.mat','.sbml','.xml'})),:)=[];

modelsToLoad={};

for i=1:length(models)
    modelsToLoad{i} = [modelFolder filesep models{i}];
end

% check if file already exists
if ~isfile([propertiesFolder filesep 'Reactions_' reconVersion '.txt']) && ~isfile([propertiesFolder filesep 'Metabolites_' reconVersion '.txt'])
    
    % print resource content
    uniqueRxns = {};
    uniqueMets = {};
    
    rxnsTmp={};
    metsTmp={};
    
    parfor i=1:length(models)
        restoreEnvironment(environment);
        changeCobraSolver(solver, 'LP', 0, -1);
        if strcmp(solver,'ibm_cplex')
            % prevent creation of log files
            changeCobraSolverParams('LP', 'logFile', 0);
        end
        try
            model = readCbModel(modelsToLoad{i});
        catch
            model=  load(modelsToLoad{i});
            model = model.model;
        end
        % collect unique metabolites and reactions resource content for refined reconstructions
        rxnsTmp{i}=model.rxns;
        mets=strrep(model.mets,'[c]','');
        mets=strrep(mets,'[e]','');
        mets=strrep(mets,'[p]','');
        metsTmp{i}=mets;
        
    end
    for i=1:length(models)
        % grab all unique reactions and metabolites
        uniqueRxns=unique(vertcat(uniqueRxns,rxnsTmp{i}));
        uniqueMets=unique(vertcat(uniqueMets,metsTmp{i}));
    end
    
    % load VMH database
    database=loadVMHDatabase;
    
    % print out the unique reactions and metabolites of the resource-only if it
    % used VMH nomenclature
    reconMetabolites=database.metabolites;
    [C,IA] = setdiff(reconMetabolites(:,1),uniqueMets);
    reconMetabolites(IA,:)=[];
    if size(reconMetabolites,1)>200
        writetable(cell2table(reconMetabolites),[propertiesFolder filesep 'Metabolites_' reconVersion],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
    end
    
    reconReactions=database.reactions;
    [C,IA] = setdiff(reconReactions(:,1),uniqueRxns);
    reconReactions(IA,:)=[];
    if size(reconMetabolites,1)>200
        writetable(cell2table(reconReactions),[propertiesFolder filesep 'Reactions_' reconVersion],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
    end
end

end
